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Abstract 

Discoveries of exoplanets orbiting evolved stars motivate critical examinations of the dy- 
namics of ./V-body systems with mass loss. Multi-planet evolved systems are particularly 
complex because of the mutual interactions between the planets. Here, we study the un- 
derlying dynamical mechanisms which can incite planetary escape in two-planet post-main 
sequence systems. Stellar mass loss alone is unlikely to be rapid and high enough to eject 
planets at typically-observed separations. However, the combination of mass loss and planet- 
planet interactions can prompt a shift from stable to chaotic regions of phase space. Conse- 
quently, when mass loss ceases, the unstable configuration may cause escape. By assuming 
a constant stellar mass loss rate, we utilize maps of dynamical stability to illustrate the 
distribution of regular and chaotic trajectories in phase space. We show that chaos can drive 
the planets to undergo close encounters, leading to the ejection of one planet. Stellar mass 
loss can trigger the transition of a planetary system from a stable to chaotic configuration, 
subsequently causing escape. We find that mass loss non-adiabatically affects planet-planet 
interaction for the most massive progenitor stars which avoid the supernova stage. For these 
cases, we present specific examples of planetary escape. 



1 Introduction 

The first confirmed extrasolar planets were found to orbit evolved stars |65[ [66] and the last 5 
years has seen a resurgence of interest in this topic due to new discoveries (e.g. [EEJ EE1 E]). 
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The majority of complimentary theoretical analyses has focused on single-planet, single-star 
evolved systems [61] E21 E7J H21 [52] or systems with a belt or disc of material [5j [61 [12]. These 
investigations crucially establish physical and analytical frameworks from which to explore more 
complex systems in greater depth. 

Just a few studies have considered multi-planet post-main-sequence (MS) evolution. [T5] per- 
formed long-term integrations of our Solar system assuming a constant Solar mass loss rate. [11] 
demonstrated that multiple equal-mass planets on coplanar, circular orbits which are marginally 
stable on the MS can become destabilised during post-MS evolution, because of the expansion 
of the Hill-stability limit due to mass loss. [58] showed that even Hill-stable pairs of planets 
can eventually become unstable, illustrating that instability occurs more readily than previously 
thought, and that planets do not need to be closely packed to become unstable many Gyr later. 
[46j considered the evolution of two planets in cataclysmic variable systems, but did not focus 
on planetary instability. 

Although these initial studies have now broached the topic of planet-planet scattering amidst 
mass loss, the detailed nature of this instability has yet to be explored. Here, we pursue this line 
of investigation. We consider a planetary system with two planets in a resonant or nonresonant 
configuration that is stable on the MS and is perturbed when mass loss takes place. We will 
show that the perturbation pushes the system into a chaotic regime through which instability 
manifests itself after the mass loss ceases. 

The consequences of this late time instability may include collisions within the system or 
escape from the system. We focus on the later possibility, due to its potential relevance to 
the purportedly vast free-floating planet population in the Milky Way |54j . which is thought to 
outnumber the Galactic population of bound planets. The existence of such substellar objects, 
known as free-floating or orphan planets (e.g. [33 [ I68 [ [13] ) . may help us understand the low-mass 
end of the initial mass function. 

The source of free-floating planets is unknown. |60j considered the rate of ejections in planet- 
planet scattering in MS systems, and found that even in the most optimistic case, the rate is 
insufficient to generate the observed poputation. Alternatively, [56J illustrated how a free-floating 
planet passing a stellar system could trigger dynamical instability and incite planetary escape, 
leading to another free-floating planet. [67] and [59J considered similar mechanisms for stellar 
flybys and star-planet flybys in the Galactic disc, respectively. These studies do not reproduce 
the required escape rate, and potentially suggest that the free-floating planet population is a 
relic of the initial formation process. The turbulent early-age birth environment of planets 
might represent the primary source, as well as the frequent and slow flybys characteristic of 
young clusters [21 [TTl [53j [34j HJ [45] . Our focus here, however, is on planetary escape in evolved 
systems, and on bounding the phase space in which this can occur. 

The paper is organized as follows. In Section 2 we present some analytic and numerical results 
for the orbital element evolution of planets under stellar mass loss without taking into account 
the mutual planetary interaction. In Section 3 we present results that show the distribution of 
chaos and order in a system consisting of two massive planets and a star with constant mass; 
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we discuss the resulting chaotic evolution and the destabilization of the system. In Section 
4 we present our results, which represent numerical simulations that combine both planetary 
interactions and stellar mass loss, and demonstrate the possible system destabilization and 
planetary escape. We present our main conclusions in Section 5. 

2 The two body problem of variable mass 

In the classical two body problem, the two bodies describe similar elliptic orbits around their 
center of mass when the motion is bounded. The motion is planar and the orbital elements of 
the two bodies (or the orbital elements of the relative motion of one body around the other) 
are constant. If, however, one of the bodies loses mass isotropically in all directions, then the 
evolution of the system is quite different. The orbital elements of the relative motion are no 
longer constant, and the evolution of the system depends on the mass loss function m = m(t), 
where m(t) is the sum of the masses of the two bodies. In this sense, the assumption of that 
mass loss by the star and mass gain by the planet are equivalent does not hold. A particular 
case, which we shall consider in this work, is a planetary system with a star and a planet with 
much smaller mass, where the star loses mass isotropically. 

The osculating elements of the relative motion are given by the differential equations [22J: 

da _ l+2ecos /+e 2 m 
dt ' 1— e 2 m 

§ = _( e + COS /)™ 

duj _ sin / rn ( ) 

dt em 

d£ (Gm) 1 1 ' 2 (1+e cos ff sin / rh 
dt [a(l-e 2 )] 3 / 2 l ~ e m> 

where a is the semimajor axis of the relative motion, e the eccentricity, u the argument of 
pericenter and / the true anomaly. One can readily check that the above equations admit the 
integral 

c 2 = Gma(l - e 2 ), (2) 

which in fact is the angular momentum integral. We set G = 1 in the remainder of the paper. 

In the following, we assume that the loss of the mass of the star is described by the Eddington- 
Jean's law 

m = —am n , (3) 

where a is, in general, a small positive constant. The exponent n defines the particular law, 
e.g. for n = we get a linear variation of mass, m = tuq — at, and for n = 1 the exponential 
one, m = moe~ at . Also, we assume that the star loses mass down to a minimum limit m m j n = 
(1 — /3)mo, where the constant (3 defines the final mass loss ratio (0 < f3 < 1). This minimum 
mass is obtained at a time t = t£, which is found through the solution of (|3|) and depends on 
mo = m(0) and the parameters a, (3 and n. 
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We introduce the normalized semimajor axis a = a/ao and the normalized mass m = m/m®. 
The evolution of the osculating elements in the interval < t < tg can be given as a function of 
the normalized mass m, which decreases monotonically in time. Then equations ([1]) are written 
as 



da _ l+2e cos f+e 



™Jl = -( e + cos /) (4) 

dm e 

df _ m 3 ~ n r-i i „ f\2 i sin/ 



m£ = -=^(l + e coB/)' + 
where 

e = ac 3 mjf- 3 . (5) 

From the above equations, and taking into account that a(0) = 1 and m(0) = 1 always, we 
conclude that two orbits with the same initial eccentricity e$ and true anomaly /o, show an 
equivalent evolution if the system is characterized by the same value of the parameter e. The 
initial value of the argument of pericenter ujq affects only the evolution of u = oj(t). 

In the above analysis, we considered a star and only one planet. If there are two, or more, 
planets, then the evolution of the osculating elements becomes more complicated, because the 
gravitational interaction between the planets affects in an important way the orbital evolution, 
as we shall see in the following sections. Only in the very special case where the masses of the 
planets could be considered as negligible would the evolution of each planet be determined by 
the system (jl]) or (jlj). 

2.1 Slow variation of mass - Evolution of the semimajor axis 

It can be proved p3] that for a slow variation of mass and small eccentricities, the action 
J = \Jma is an adiabatic invariant, which implies that 

m(t)a{t) rj constant. (6) 

Thus, from the angular momentum integral ([2]), we find that the eccentricity e remains constant. 
From the above we see that the secular variation of the semimajor axis is given by 

m{t) 

which implies that the semimajor axis increases monotonically, since the mass decreases. Par- 
ticularly in the case of linear mass loss, m = mo — at, we obtain the adiabatic estimation 
a-adiab = o,o(l — at/m®) -1 given by Veras et al. (2011), while at the minimum mass limit the 
semimajor axis takes the value 

a £ = a /(l - /3). 
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2.2 The evolution of eccentric anomaly 



The variation of the eccentric anomaly E in time, when the star loses mass isotropically, is 
necessary for estimating the time evolution of the eccentricity and the argument of pericen- 
ter, as we will show in the following paragraphs. The eccentric anomaly obeys the equation 
(Hadjidemetriou, 1966) 

dE 1 — e cos Edf s'mE de 

~dT = (1 - e 2)i/2 dt ~ T^dt' () 

where the derivatives of the true anomaly and the eccentricity are given by equations ([I]). In 
order to find the secular change of E, we ignore in ([8]) all the periodic terms and then, up to 
first order terms in eccentricity, we obtain 



dE m 2 
~dt ~ ~cF' 



(9) 

or, by considering the normalization used in equations 



m— = , (10) 

dm e 

From the above equations we obtain that dE/dt > always holds true. The evolution of E as 
a function of m depends essentially only on the system's parameters e and n. For the case of 
linear mass loss (n = 0) we obtain 

^ „ 1 — 
E(fh) = Eo^^, 

where Eq = E(0), or 

E(t) = E + ^(3m (m -at) + a 2 t 2 ). (11) 



2.3 Analytic estimate of the evolution of the eccentricity and the argument 
of pericenter 

Hadjidemetriou (1966) has given a series solution of equations (pQ) for the eccentricity and the 
argument of pericenter as a function of the eccentric anomaly E, 

e = e + ^2e J fj(E;e ) 
i=i 

w = uj + ^V#j(-E;e ), 
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Figure 1: Evolution of eccentricity for the initial conditions ao = 1, eo = 0.05, /o = and linear 
mass loss with tuq = 1 and a = 0.001. The dashed line corresponds to the numerical solution of 
equations ([1]). The solid bold line is the analytic solution (|12p . where E have been substituted 
by eq. (Illj) . The estimate (I13p of the maximum eccentricity is also shown by the thin solid line. 

which converge for e ^ ^2^3/2 < 1- Assuming e as a small parameter we can write, for small 
eccentricities, an approximate solution for e = e(E) up to 0(e 3 ) as 

e 2 

e sa eo + esin^H sin 2 E + e 2 E(3 — n) sin E 

2e 

+ e z E 2 (3-nfsinE + 0(e i+1 E i ), i > 3. (12) 

The evolution of the eccentricity in time is given by considering the particular solution E = E(t) 
of Eq. 0. Apart from the constant term eo, the first two terms in (|12p are periodic terms of 
period 2tt (one planetary revolution). The last two terms indicate also the same periodicity but 
with a secular variation of the amplitude of the periodic oscillations. Thus, an estimate of the 
maximum value reached by the eccentricity along the evolution can be approximated by the 
formula 

A e (E) = A + e 2 \3-n\E + e 3 (3-n) 2 E 2 + 0(e t+1 E i ), i > 3. (13) 

The remaining terms 0(e l+1 E l ) are small for the first planetary revolutions, but they become 
important as the number of revolutions increases. In Fig. (prj we show an example of the 
evolution of eccentricity. We see that the analytic solution (|12p coincides very well with the real 
(numerical) solution but only for a relatively short time interval. 

Working in the same way as with eccentricity, we can write an approximate solution for the 
argument of pericenter 

9 n — 3„ q o (n — 3) 2 „ 

oj w ojq + e cos E + e 2 E cos E - e 6 E 2 - ^cos£ 

eo e 
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Figure 2: The evolution of the eccentricity (left panels) and the argument of pericenter (right 
panels) as the star loses mass with a linear rate, m(t) = mo — at. Different color lines correspond 
to different initial values eo indicated in panel (a). 



+ Wa + Oie^E'), i>3), (14) 

where 



Wa = e^E-e* n{n - 3) E* 
2 2 



is a secular component of the evolution. Note that w a is zero in the case of a linear mass loss 
(n = 0). 



2.4 Numerically calculated orbital evolution for one planet 

The analytic estimates found above are sufficiently accurate only for time intervals where the 
mass loss does not exceed 20-30% of the total mass. However, the numerical solution show that 
the eccentricity e and the argument of pericenter uj continue to show oscillations with increasing 
amplitude for still larger mass loss, until their evolution enters a "runaway regime" where the 
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Figure 3: The evolution of the eccentricity for different exponents n in the Eddington-Jean's 
law given by Eq. ([3]). For each case we give the evolution for e = 0.001, 0.01 and 0.1, as it is 
indicated in panel (d). The initial conditions are eo = 0.1, /o = ujq = 0. 



eccentricity will eventually increase monotonically. The passage of the evolution to the "runaway 
regime" has been studied in detail by Veras et al. (2011) for the case of a linear mass loss. In 
the following we present some results obtained by numerical integrations. 

In Fig. [2] (panels a,b and c) we show some typical examples of the evolution of the eccentricity 
e, where the horizontal axis represents the mass loss ratio (3 = 1 — m, which can be mapped 
to the time t through the particular mass loss law of eq. ([3]). Here we consider the linear law 
(n = 0). We observe that when the orbits enter the "runaway regime", we soon get e > 1 and 
the planet escapes. Such an escape is obtained when the mass loss becomes sufficiently large, 
equal to a critical value ft*. This value depends generally on the initial conditions of the system, 
but when we start with small eccentricity values eo, /3* depends mainly on the parameter e. We 
observe that for e = 0.1 the orbit enters the runaway regime from the beginning. In the right 
panels d, e and f of Fig. [2] we present the evolution of the argument of pericenter uj of the same 
orbits as in the left panels. Again we find that for small values of e, oj librates with an increasing 
amplitude, as it is indicated by the approximation of eq. f|14|) . However, in the case of e = 0.1 
(panel d), oj increases monotonically. In all cases oj completes one revolution at most. 

By considering n / 0, i.e. non-linear mass loss, we find that the evolution of the orbital 
elements is qualitatively similar to the linear case even though the runaway regime appears in 
general for larger values of mass loss. Some examples of the evolution of the eccentricity are 
shown in Fig. [3] for the cases n = 1/2, n = 1 (exponential decay), n = 2 and n = 3. For the 
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case n = 2 (panel c) and for e = 0.001, the runaway behaviour appears only when the star mass 
becomes very small. An exception is the case n = 3 (panel d) where the eccentricity continues 
to show small oscillations as m — > and the runaway regime does not exist (at least for any 
e < 0.1). From the series solution Eq. (|12p we obtain that for n = 3 the amplitude of the 
oscillations is constant and the eccentricity is a periodic function of the eccentric anomaly (see 
also Hadjidemetriou, 1966). 

So far we have shown that the planetary eccentricity may take large values either due to 
its large amplitude oscillations or due to the entrance in the runaway regime. Let us assume 
that the star mass loss starts at t = 0, when e = eo, and continues up to a particular mass 
loss ratio value /?, which is reached after time t = t^. Then, at the end of the star mass loss 
event, the planetary orbit has reached a new eccentricity value et = e(te). In the contour maps 
of Fig. U] we present the value of in the parameter domain < ft < 1 and 0.001 < e < 1. 
The cases of linear (left panel) and exponential (right panel) rates of mass loss are given. The 
tabulated eccentricity value is the average eccentricity obtained at t = from 100 orbits with 
initial eccentricity and true anomaly values that were randomly selected in the intervals (0,0.1) 
and (0,27r), respectively. We see that for small values of the parameter e (e.g. e < 10 -3 ), a 
significant increment of the eccentricity is obtained only after a significant mass loss (e.g. more 
than 80% in the linear case and more than 90% in the exponential case). As e takes larger 
values, a significant eccentricity increment is observed for smaller amounts of lost mass. The 
escape regime (e > 1) is indicated by the gray shaded region and occupies the region of large 
values of f3 and e, as was expected. By computing the same maps for n = 2 and n = 3, we 
obtain similar plots. However, these curves, which correspond to the same eccentricity levels, 
are obtained for larger (3 values and yield smaller escape regimes. 

2.5 Resonant evolution 

As previously mentioned, when there are at least two planets in the system, the evolution 
does not depend on the mass loss of the star only; the gravitational interaction between the 
planets plays an important role. In order to study this effect, and reveal the extent to which 
the gravitational interaction influences the evolution, we consider here a planetary system with 
two planets of negligible mass amidst mass loss. In this case, no gravitational interaction exists 
between the planets. 

Consider a system with two massless planets p, i = 1,2, with semimajor axes and 
angular momenta c%. The evolution of each planet depends on the value of the parameter 
e, namely = ac^mg" 3 . Because is different for each planet, each planet then follows a 
different evolution. If Pi is the outer planet (02 > ai), then we obtain €2 > e\. The mean 
motion resonance is defined by the mean motion ratio n\/n2, where ni = m 1 / 2 a i 3//2 , which is 
given by the relation 




(15) 
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Figure 4: The mean eccentricity values obtained when the star mass loss ceases at a value 
m = mo(l — /3) and for e in the range (0.001,1.0). The mean eccentricity is obtained after 
running 100 orbits with different initial conditions, < eo < 0.1, < /o < 2n and u>q = 0. The 
gray region indicates the escape regime, where e > 1. Left and right panels correspond to the 
linear (n = 0) and the exponential (n = 1) rate of mass loss, respectively. 



Thus two resonat planets with orbits of low eccentricities preserve their mean motion resonance 
^ 2^2, p } q integers, under the star mass loss and before entering the runaway regime, 
where the eccentricities take large values. Also, the preservation of the resonance and the 
libration of the argument of pericenter, mentioned in section \2A\ indicate that the resonant 
angles Aw = uj 2 — uj\ and 9^ = pXi — (p + q)X 2 + qoJi, i = 1,2, should also librate, until a 
significant stellar mass loss is reached. 

In Fig. [5] we present the evolution of the mean motion ratio and of the resonant angles for 
the 2/1 (left panel) and 3/1 (right panel) resonance and for e\ = 10~ 3 , 10 -2 and 10 _1 . For the 
outer planet it is e 2 « 2ei and e 2 ~ 3ei for the 2/1 and 3/1 resonance, respectively. For both 
resonances we see that the system evolves in a similar way. For small mass loss parameters e,, 
n\jn 2 is almost constant before a significant amount of stellar mass loss. First the outer planet 
enters the runaway regime and then a rapid increment of n\ln 2 is observed. As we consider 
larger values of e, the system leaves the resonance at a smaller mass loss ratio f3 = 1 — fh. As far 
as the system remains close to the resonance, the resonant angles show small oscillations around 
their initial value. Afterwards, the amplitude of the libration increases and the evolution of the 
resonant angles might end up in circulation. For large values of e (e.g. for e = 0.1), we obtain a 
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Figure 5: The evolution of the mean motion ratio and the resonant angles Au and 9 = 9\ as 
the star loses mass at a linear rate, m{t) = tjiq — at. (a) starting from the 2/1 resonance and 
(b) starting from the 3/1 resonance. The initial conditions are eo = 0.1, /o = ujq = 0. 



slow monotonic variation of the angles. Particularly, Aw decreases while 9\ increases. We stop 
the evolution when the outer planet escapes, i.e., when &i > 1, where the mean motion is no 
longer defined. 

3 Chaos and ejection in two-planet systems 

Here we study the evolution of a two-planet system where the mutual planetary gravitational 
interaction is included but the star does not lose mass. Thus we establish a dynamical basis 
with which to compare the mass loss case in the next section. 

3.1 Model and methods 

The study of the dynamics of a system consisting of a star of mass mo and two planets Pi, 
i = 1,2 of masses uii can be studied by using the model of the general three body problem 
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(GTBP) 

f i = G Z) ! T* r V' * = 1,2), (16) 

i=o 

where indicates the position of the planet and = — rj. We restrict our study to the 
planar case, rj = (xi,yi) where the inertial frame Oxy is centered at the barycenter. Also we 
normalize the units by setting G = 1 and tuq + mi + mi = 1. In the following, we should take 
always the P\ to be initially the inner planet and P2 the outer one in the sense that a\ < 02, 
where is the semimajor axis of Pj. 

Although the system appears having four degrees of freedom, we can use the angular mo- 
mentum integral and a rotating frame and reduce the system to three degrees of freedom (Had- 
jidemetriou, 1975). In such a system, regular planetary orbits correspond to quasiperiodic 
trajectories, which twist on invariant tori in phase space according to the KAM theorem. How- 
ever the mutual planetary interaction destroys the integrability of the system and chaotic orbits 
coexist beside the regular ones. 




1 2 3 4 5 

log-io(t) 



Figure 6: The LCN evolution for a regular orbit (Oi), chaotic orbits (O2, O3) and a chaotic 
orbit (O4) that shows a close encounter after 3.5 Ky (we note that 6.28 time units correspond 
to about 1 year). The initial conditions are a\ = 1, 02 = 1.74, uj\ = 0, uji = 180°, Mi = M2 = 0, 
ei = 0.1 and e2 = 0.05, 0.1, 0.25 and 0.26 for the orbits O1-O4, respectively. 

Many numerical tools for the detection of chaos have been proposed. For example, in plane- 
tary dynamics the chaos indicators MEGNO, RLI and FLI have been used (Gozdziewski, 2005; 
Sandor et al, 2007; Voyatzis, 2008, respectively). In the present work we compute the maximal 
Lyapunov characteristic number (LCN), which is the classical measure for the average expo- 
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nential divergence of nearby orbits. An example of the evolution of LCN for some different 
trajectories of the GTBP is shown in Fig. [6] for 10 6 time units (or ~ 150 ky). A regular evolu- 
tion is described by an LCN evolution that tends to zero, as e.g. the orbit 0\. For chaotic orbits 
LCN tends to a positive value, as e.g. in the orbits O2 and O3. This value is used to form dy- 
namical stability maps (DS maps). In these we define plane grids of initial conditionspresented 
with a color scale representing the LCN value after a particular integration time tmadj- For the 
numerical integration of Eq. (|16p we use the Bulirsch-Stoer integrator. When planetary close 
encounters occur, the integration may break, in the sense that the integration step becomes very 
small in order for the method to preserve the requested integration accuracy (e.g. the orbit O4). 
Such cases always correspond to strongly chaotic orbits and in the DS maps are presented by 
the value LCN = 1 (light-coloured regions). 



10- 1 s 




6 0.4 0.8 1.2 1.6 2 

t (years) 



Figure 7: An example of a planetary evolution (1/1 resonant) that shows a sequence of close 
encounters a) the planetary distance ryi b) the evolution of the energy error along the numerical 
integration with and without regularization (solid and dashed curves, respectively). 

In cases where we need to follow the evolution after a close encounter and without loss 
of accuracy, we apply a regularization of the collision singularities. In this paper we apply a 
Levi-Civita transformation when the gravitational interaction between any two bodies of the 
system becomes very strong (Marchal, 1990; Aarseth, 2003). In Fig. Owe present an example 
of evolution with a sequence of planetary close encounters. The two planets are placed inside 
the Hill radius and evolve in a satellite configuration. In panel (a) we present the planetary 
distance r\i- We integrate the system using two methods. In the first we integrate equations 

An application of LCN maps of dynamical stability to planetary dynamics is given in Hadjidemetriou and 
Voyatzis (2011a). 
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Figure 8: DS maps in the eccentricity plane for n\/n 2 ~ 2.29 (a\ = 1.0, 02 = 1.74) (panels a,b) 
and ni/ri2 ~ 3.29 (ai = 1.0, a 2 = 2.21) (panels c,d) . The top panels correspond to m\ = mj, 
m 2 = 0.2m j and the bottom ones to mi = mj, m 2 = 2m j. For all orbits we set as initial 
conditions w a = 0°, co 2 = 180°, M x = M 2 = 0°. The dashed curve indicates the collision line. 



(|16|) using the BS integrator with accuracy 10~ 14 . When the step becomes very small, instead 
of stopping the integration, we reduce the requested accuracy. In the second method we use 
the same integrator as above but we switch on the regularized equations during planetary close 
encounters. In panel (b) of Fig [7] we present the error in energy along the trajectory in both 
cases. We can observe that without regularization the error increases after each close encounter, 
while the integration of the regularized equations preserves the accuracy at a very good level. 

3.2 DS maps 

We restrict our numerical simulations to planetary systems with the inner planet having mass 
mi = 0.001 (approximately Jupiter's mass, mj) and semimajor axis a\ = 1 AU. The outer 
planet is either lighter (m 2 = 0.2mj) or heavier [m 2 = 2mj) than the inner one and with 
semimajor axis in the range 1.5 < a 2 < 2.5 AU. 
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3.2.1 Non-resonant motion 

For non-resonant motion, the main source of chaos is planetary close encounters, which can 
happen when the initially elliptic orbits of the two planets intersect. In Fig. [8] we present DS 
maps for non resonant motion. The maps are defined by grids of initial conditions in the plane 
of eccentricities. The other orbital elements are fixed and are given in the caption. In panels 

(a) and (b) the orbits correspond to an initial mean motion ratio n\jni = 2.29. In case (a), 
where the outer planet (P2) is more massive than the inner one (-Pi), we observe that regular 
orbits exists only in a well-defined region at small eccentricity values. Some small stable regions 
appear for e2 up to 0.3 but only for e\ < 0.1. A similar picture is obtained also for the case 

(b) , where the inner planet is lighter than the outer one. Now the stable orbits are confined in 
a region of small eccentricities. Above the line a\{\ + e\) = 02(1 — &%) (collision line), which is 
represented by the dashed line in the dynamical maps, the planetary orbits intersect and strong 
chaos appears. We obtain the same dynamical picture in panels (c) and (d), where n\/ri2 = 3.29. 
In this case the outer planet starts with a larger semimajor axis than in the previous case. The 
planetary interactions are weaker, the stable region covers a larger domain and the collision line 
appears at higher eccentricity values. 

3.2.2 Resonant motion 

In Fig. [9] we present DS maps for n\/n2 = 2.0 (panels a and b) and ni/712 = 3.0 (panels c and 
d). In these cases we obtain a complex structure in phase space where regions of stable and 
unstable motion are intermingled (Michtchenko et al, 2008a, b; Hadjidemetriou and Voyatzis, 
2009). Since resonances may offer phase protection mechanisms, stable long term evolution 
might appear even for intersecting orbits. We observe in the DS maps that we can have regular 
evolution when the orbit of the inner planet is very eccentric. However when the eccentricity of 
the outer planet becomes larger than about 0.2 we obtain a wide chaotic sea where trajectories 
are strongly chaotic. For the 3/1 resonant case with m-i > rri2 we observe a zone of regular 
orbits for 0.25 < e2 < 0.35. The existence of such zones and, generally, the distribution of 
chaotic and regular regimes in phase space is significantly affected by the families of periodic 
orbits (symmetric or asymmetric) that, generally, exist in the resonant regions (see e.g. Voyatzis, 
2008). 

3.2.3 Exact resonances - periodic orbits 

Resonant motion is related with the existence of monoparametric families of periodic orbits in 
phase space. Such periodic orbits are also called exact resonances (Beauge et al, 2003), which can 
be stable or unstable and contribute significantly to the topology of phase space and the existence 
of order and chaos. In particular, stable periodic orbits are surrounded by regions of regular 
motion where the resonant angles Acj, 6\ and 82 librate. Outside these regions stable motion 
can be also obtained, where some of the resonant angles circulate (Michtchenko et al, 2008a, b; 
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Figure 9: DS maps on the eccentricity plane for the 2/1 and 3/1 resonances. The top panels 
correspond to m\ = mj, m 2 = 0.2mj and the bottom ones to m\ = mj, m 2 = 2m j. The initial 
semimajor axes are a\ = 1.0 and 02 = 1.59 and 2.08 for the 2/1 and 3/1 resonances, respectively. 
For all orbits we start with uj x = 0°, uj 2 = 180°, M x = M 2 = 0°. The dashed line indicates the 
planetary collision line. 
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Figure 10: DS maps around a 2/1 resonant stable periodic orbit on the plane of resonant angles 
Aw). The periodic orbit is located at (0,0) with (approximate) initial conditions a) a\ = 1, 
a 2 = 1.582, e\ = 0.3, e 2 = 0.063 b) a\ = 1, 02 = 1.580, e\ = 0.5, e 2 = 0.166 and c) a\ = 1, 
a 2 = 1.580, ei = 0.7, e 2 = 0.318. For all cases it is wi = w 2 = 0, Mi = 0, M 2 = 180° and 
mi = mj, m 2 = 2m j. The cross in the center indicates the position of periodic orbit. 

Voyatzis, 2008). In Fig. [10] we present maps of dynamical stability which show the distribution 
of regular motion around 2/1 resonant periodic orbits, which belong to the symmetric family of 
periodic orbits for mi = mj and m 2 = 2m j with configuration (9\,Auj) = (0,0) (family Si in 
Voyatzis et al, 2009). The maps are defined on the plane of resonant angles (0i, Aw) and we see 
that the regular region, which surrounds the periodic orbit at (0, 0), shrinks as the eccentricities 
increase. 

Resonant families of stable periodic orbits are of particular importance, since they form 
paths in phase space for a migrating planetary system (Lee and Peale, 2002; Ferraz-Mello et 
al, 2003, Hadjidemetriou and Voyatzis, 2010, 2011b). So, regions around stable periodic orbits 
are candidates for containing a planetary system after a migrating process and a capture in 
resonance. 

3.3 Escape of planets from planetary systems in a chaotic region 

In this section we study the evolution of a planetary system under the gravitational interaction 
between the two planets, when the system is initially located in a chaotic region. The mechanism 
that transfers an initially stable system to a chaotic region will be studied in the next section. 

Since the system is of more than two degrees of freedom, invariant tori are not boundaries 
for the chaotic regimes in phase space and fast or slow diffusion of chaotic motion is possible (see 
e.g. Lichtenberg and Lieberman, 1983). Also, if we take into account that in the GTBP zero 
velocity curves do not exist to bound the planetary orbits (Marchal and Bozis, 1982), then we 
may conjecture that chaotic motion leads to the escape of a planet after a long term evolution. 
This has been established after extensive numerical simulations of the GTBP (see Valtonen and 
Karttunen 2006, and references therein). In the following we discuss the particular case of a 
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Figure 11: The mutual distance of planets P\ and P2 along a strongly chaotic orbit with initial 
conditions as those in the DS map of figure [9^1 and for e\ = 0.05 and ei = 0.3. The planet P2 
escapes after about 30 Ky. 

planar two-planet system. 

In the dynamical maps presented in the previous subsections, we observed wide chaotic re- 
gions with strongly chaotic motion. Many numerical integrations of orbits with initial conditions 
inside these regions showed that, during the evolution, close encounters between planets occur 
and the less massive planet is scattered to orbits with larger eccentricities. The continued irreg- 
ular evolution results in a sequence of close encounters that ejects the planet to large distances 
and finally to escape. We classify the evolution as escape when a planet that starts from a 
distance of order 1 AU, moves to a distance larger that 1000 AU with eccentricity e > 1. In 
these numerical integrations we set the accuracy in the BS integrator to 14 digits and we use 
the regularized equations during close encounters. 

In Fig. [11] we present the evolution of the planetary distance for an orbit starting from the 
strongly chaotic (light colored) region appearing in the map of dynamical stability of Fig. [9^ at 
(ei,e2) = (0.05,0.3). The first close encounter appears after about 8 Ky. A sequence of close 
encounters follows and after 20 Ky scattering of the planet P2 to large distances is observed. 
Finally the planet escapes at about 30 Ky. In Fig. [T2k we present the evolution of the semimajor 
axis and the eccentricity. After the first close encounter the eccentricity and the semimajor axis 
of P2 show a jump to higher values and their variation in time becomes very irregular. The orbit 
of the heavier planet Pi shows also the same irregularity, but the variation of its orbital elements 
is quite smaller. Finally we get e2 > 1 and for the remaining planet e\ = 0.13, a\ = 0.89. 

A second example of chaotic evolution is presented in Fig. [12b. Now we start from the 
point (ei,e2) = (0.05,0.26) which is closer to the regular (dark colored) region of the DS map 
of Fig. [9^,. Up to 120 Ky the evolution of the semimajor axis and eccentricity of both planets 
seems vary regular. Afterwards chaotic evolution appears and a sequence of close encounters 
takes place, similar to the previous case. Finally planet P2 escapes and the remaining planet Pi 
moves on an elliptic orbit with e\ = 0.12 and a\ = 0.89. By considering initial conditions that 
correspond to chaotic motion but are closer to the regime of regular motion, we obtain similar 
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Figure 12: The evolution of semimajor axis and eccentricity of the planets P\ and P2 along 
a chaotic orbit with initial conditions same to those in the DS map of figure [9^t and initial 
eccentricities (a) e\ = 0.05, e2 = 0.3 (b) e± = 0.05, &2 = 0.26. 

planetary destabilization, which, however, occurs after a longer time span. For example, by 
considering the same initial conditions as above but for e2 = 0.25 and e2 = 0.23 we found that 
escape appears at 0.9My and lOMy, respectively. However we should remark that the escape 
time is very sensitive to the numerical integration accuracy. 

Chaos is a necessary but not sufficient condition for planetary ejection over reasonable time 
spans (e.g. MS lifetimes). If we consider the initial conditions in a narrow chaotic zone in the 
maps of dynamical stability and far from the wide chaotic sea, then chaos may remain bounded 
throughout long-term evolution. In other words, consider a system with initial conditions in the 
chaotic region comprising small eccentricities (ej < 0.15) of the DS map of Figj9k. Although 
this system evolves apparently irregularly, up until lOOMy the system is stable in the sense that 
no ejection or collision occurs. The existence of bounded chaos for orbits which are Hill stable 
is also noted by [20j. 

4 Evolution under mass loss from the star and gravitational 
interaction between the planets 

Here we present our main results. We analyze how mass loss could couple with the planets' 
mutual interaction to cause one of them to escape. In the following study we start with a 
two-planet system, which is trapped in a stable (resonant or non-resonant) configuration. Then 
we assume that the star begins to lose mass in an isotropic way, but in a manner where the 
percentage of the mass lost and the rate of mass loss are not large enough to cause the direct 
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Figure 13: The evolution of the planetary semimajor axes and the eccentricities for the periodic 
orbit of Fig. PTOk when it is affected by the star mass loss with (3 = 75% and a = 1CP 3 . In panel 
(a) the first stage of the evolution is presented where the vertical dashed line indicates the time 
ti where the mass loss ceases. In panel (b) the continued evolution is presented until the escape 
of the planet Pi . 



ejection of the planet. However, as we showed in Section 2, even in these cases the orbital 
elements of the planetary orbits are affected. So the system migrates in phase space and after 
the end of the mass loss process it can be found inside a chaotic region. From that point on, the 
chaotic evolution of the system suffers from close encounters, allowing for the possible ejection 
of a planet (see section I3.3|) . 

We consider a linear rate of mass loss, m = mo — at, which takes place up to time tg = f3mo/a 
where the mass of star is reduced to m = (1 — /3)mo (see section 2). In all the figures showing 
the evolution of the system, the inner planet has the subscript "1" and is represented by a blue 
line and the outer planet has the subscript "2" and is represented by a red line. 

4.1 A fiducial example 

Here, the initial total mass of the system is set to one solar mass, namely niQ ~ 1. As an 
example we consider the stable 2/1 resonant periodic orbit given in Fig. [10a and we apply a 
star mass loss /3 = 75% and rate a = 10~ 3 . The evolution of the planetary semimajor axes 
and eccentricities for the first years of the mass loss is shown in Fig. [T3"a . The semimajor 
axes increase for both planets approximately according to the adiabatic estimate given by Eq. 
([7]), but the ratio of mean motions n\/n2, represented by the horizontal dashed line, is almost 
constant, as it is suggested by Eq. (Tl5l) . After the end of the mass loss process, which ceases 
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Figure 14: The eccentricity evolution of an orbit located near the periodic orbit of the DS map 
of Fig. [TUh . The initial conditions are those of the periodic orbit except that we set UJ2 = 40° 
and the mass loss parameters are a = 10 -4 and j3 = 50%. 

at 120 years, the eccentricities have increased slightly but the mass of the star has decreased 
significantly. Such new conditions corresponds now to a chaotic orbit as it becomes evident in 
Fig. 113b . We observe that the eccentricities incur well-bounded but weakly chaotic oscillations 
up to 70 Ky. Then the system enters a strongly chaotic region and the eccentricities show large 
and irregular variations up to about 290 Ky. In this interval, the variation of the semimajor axis 
of the planet P2, which is the heavier one, is significant. After this time interval the semimajor 
axis and the eccentricity of the planet Pi increase rapidly and finally the planet is ejected. It is 
obvious that if we ignore the gravitational interactions between the planets, the orbits of both 
planets would be constant ellipses with orbital elements attained at t = ti. 

The above is a typical example of destabilisation after stellar mass loss. The same results 
are obtained for the periodic orbits in panels (b) and (c) of Fig. [TUl As we mentioned in 
section 3, stable periodic orbits are in the centers of islands of stability. If we consider initial 
conditions far from the periodic orbit and near a chaotic region, destabilization can happen for 
smaller values of the parameters a and (3. An example is shown in Fig. Q3] where we use the 
same initial conditions except that 102= 40°, which means that the orbit is located at the point 
(6>i, Acj)=(-80°,40°) of the DS map of Fig. [TDk. Now the system is destabilized for a = 10~ 4 
and p = 50%. 

4.2 Maps of destabilization 

In section [3721 we presented DS maps on the plane of eccentricities, which show the distribution 
of regular and chaotic orbits when the mass of the star is constant. According to the mechanism 
described above, the regular orbits can become chaotic if stellar mass loss takes place. In order 
to make a more extensive study of the possibility of planetary destabilization, we consider the 
initial conditions of the regular orbits in the DS maps of figures [8] and [9] and we follow the 
evolution of the system by considering star mass loss with a = 10 -4 and [3 = 75%. For each 
orbit we compute the LCN and thus we determine which of the orbits become chaotic after the 
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Figure 15: Maps which show the orbits that remain stable or destabilized after the mass loss of 
the star (a = 10 -4 , /3 = 75%). The maps correspond to the DS maps of figures [8] and EJ as it is 
indicated in each one of the maps. The ratio of mean motion r is also indicated while the index 
(i) or (e) indicates that the less masive planet is the inner or the external one, respectively. The 
black areas indicate the orbits that continue to evolve regularly after the effect of the star mass 
loss. The gray area indicates the stable orbits that became chaotic. The percentage of the orbits 
that became chaotic is indicated by the index d. 



stellar mass loss. We assume that if LCN is greater than 5 x 10 -5 , then the evolution is chaotic. 
The results are presented in Fig. [15j The black shaded areas denote the orbits which remain 
regular after the mass loss process. The gray shaded area denotes the regions where the initially 
regular orbits become chaotic. 

For each map we indicate the percentage d of the orbits that become chaotic. In the resonant 
case with r = 2 and the nonresonant one with r = 2.29, we find that a large fraction of orbits 
is destabilised. In contrast, at the 3/1 resonance and for the nonresonant case r = 3.29, the 
majority of the orbits remain regular, especially when the outer planet is the less massive one. 
We note that in the first cases (r = 2 and r = 2.29) the planets are closer to each other and 
their gravitational interaction is stronger. Thus, we may claim that planetary orbits are more 
affected by stellar mass loss as the planetary interactions becomes stronger. This important 
claim is not directly assessed in [TT] nor [58], as the simulations in those studies are not set up 
to address this issue. Further, in those studies, the planets are of the same mass and the outer 
planet is the one that escapes. In our integrations the planet which escapes is always the lighter 
one (see also [T5]), 
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Table 1: Thermally Pulsing Asymptotic Giant Branch Properties 



Progenitor Mass (3 t^/(ky) a 



8M Q 


77.8% 


78.0 


1.27 x 10- 


-5 


7M Q 


78.5% 


81.1 


1.08 x 10- 


-5 


6M 


78.8% 


89.1 


8.44 x 10- 


-6 


5M 


78.0% 


98.7 


6.26 x 10- 


-6 


4M 


76.3% 


108.2 


4.50 x 10- 


-6 


3M 


73.3% 


184.6 


1.90 x 10- 


-6 


2M Q 


65.5% 


272.1 


7.64 x 10- 


-7 



4.3 Realistic cases 

Although this study is focused on the dynamical properties of the general three-body problem 
with mass loss, we can relate the results to real systems. Doing so helps us determine in what 
cases will significant mass loss affect planet-planet scattering in a non-adiabatic manner. 

In order to make this relation, we use the SSE (Single Star Evolution) code [29] to cre- 
ate stellar evolutionary tracks. For Solar metalicity stars, and a Reimers mass loss coeffi- 
cient of 0.5, we trace the time evolution of stars with the following progenitor masses: 1M , 
2M ,3M ,4M ,5M ,6M ,7M and 8M . We find that for all cases except the 1M Q case, the 
vast majority of the mass is lost on the Thermally Pulsing Asymptotic Giant Branch (TPAGB) 
phase. This phase, which is typically one of the shortest in duration, also typically features the 
greatest mass loss rates. This combination is well-suited for observing the type of phenomenon 
seen in the fiducial simulations. Therefore, we neglect the other phases of evolution and assume 
that all the most lost occurs during the TPAGB. 

We summarize the output of SSE, which includes j3 and ti (for the TPAGB phase only) in 
Table 1. Then we remove the assumption that tuq ~ 1 and compute a and record that value 
as well in the last column. The table suggests that adopting a ~ 10 -5 may reflect real systems 
with progenitor stellar masses of 6M — 8M . 

In Fig. [T6l we show a typical example of the evolution of such a multi-planet system with 
a star of mass 8M and two planets with m\ = 3Mj and mi = 0.6Mj. The planets are 
initially located in the 2/1 resonance with initial eccentricities e\ = 0.1 and ei = 0.25. The 
numerical integration shows long term stability if the star does not lose mass. In particular, the 
eccentricities oscillate regularly about their initial values. If we introduce a mass loss process 
with a rate a = 10~ 5 and f3 = 60% then we obtain the evolution given in Fig. [TEh . In the time 
interval, t <ti the eccentricities show an adiabatic variation besides their fast oscillations. When 
the mass loss ceases, the system continues to evolve regularly with the eccentricities oscillating 
(with fast and slow components) about the values 0.13 and 0.19 for the inner and the outer 
planet, respectively. 
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Figure 16: The eccentricity evolution for a 2/1 resonant planetary system with mo = 8M© and 
mi = 3Mj and 7772 = O.QMj. The star loses mass linearly with rate a = 10~ 5 and (3 = 60% and 
/3 = 75% (panel (a) and (b), respectively). The dashed vertical line indicates the te value, where 
the mass loss ends. The initial conditions are a± = 10.0 and a 2 = 15.97, e\ = 0.1, e 2 = 0.25, 
Ul = 0°, uj 2 = 180°, Mi = M 2 = 0°. 



The evolution becomes very different if we consider a larger amount of mass loss. Particularly 
for (3 = 75% we get the evolution shown in Fig. 116b . Now, at t = ti the planetary system 
has entered a chaotic region and its evolution for t > t# is evidently chaotic. After 4.1 My 
close planetary encounters destabilize the system and the mutual gravitational attraction of the 
planets forces the outer, lighter, planet to be ejected from the system after about 5 My. 

Evolution in resonance is not a necessary condition for destabilization. Figure [T7] shows 
the evolution of a non-resonant system (771/71-2 — 2.415) which consist of a 4M Q star and two 
planets with mi = 0.6Afj, 7712 = 3Afj, a\ = 10.0, a 2 = 18.0, e\ = 0.1 and e 2 = 0.2. Note that 
in this case the inner planet is the smaller one. Now we consider realistic values of (3 = 76% 
and a = 5 x 10 -6 (see Table 1). Under the stellar mass loss the semimajor axes increase to 
ai = 40 AU and a 2 = 75 AU. 

The system shows some strong irregularities in the evolution of eccentricities and the orbit 
of the inner planet becomes very eccentric. Particularly, in the time interval 5 < t < 8 My the 
eccentricity of the outer planet oscillates about its initial value but the orbit of the inner planet 
shows eccentricity values in the interval 0.6 < e\ < 0.8. For t > 8 My close planetary approaches 
force the inner and lighter planet to become (temporarily) the outer one in the system. The 
evolution becomes very irregular and, finally, the lighter planet is ejected out to a distance larger 
than 1000 AU. We have also performed integrations (up to 100 My) indicating that, without 
mass loss, this system is stable and regular. Therefore, mass loss is the trigger for the instability 
in this 4M case. 

In general the numerical simulations show that the position of the system in phase space 
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determines essentially its robustness to perturbations caused by stellar mass loss. Also, the 
planetary destabilization seems to be sensitive to the parameters a and j3. Concerning the case 
of systems of more than two planets and also taking into account the results of [15] and 
we may claim that these systems may be more sensitive to small perturbations than 2-planet 
systems. As soon as more planets are included in a system, the regions of stability in phase 
space is reduced and instabilities are more likely to occur after a potential decrease of the star's 
mass. 

4.4 Link with Observed Exosystems 

Our study might help explain discrepancies in the observed populations of 2-planet systems 
on the MS with those beyond the MS. To date, no planets orbiting white dwarfs have been 
confirmed, although white dwarfs are present as distant companions in known planetary systems 
(G186: @8l SO] ; HD 27442: [3 El SU HI] ; HD 147513: [371 EH E] ) ■ Alternatively, over 30 planets 
are known to orbit giant stars (see Table 6 of [19]) and well over 600 are confirmed to orbit MS 
start® 

Among the mechanisms invoked to explain the general lack of known post-MS planets are 
different modes of formation around giant stars [10] and direct engulfment into the star itself 
due to star-planet tides (e.g. [HI [62j H31 [301 S2J Ej). We are proposing a third mechanism: 
scattering in two-planet systems that is induced by mass loss. The planets in these systems 
are far away enough from their parent star to remain unaffected by tidal effects. A potential 
fourth mechanism, albeit one which will be complex to model, is two-planet scattering where the 
planets are close enough to the star to be affected by both mass loss and tides. This represents 
a potentially relevant extension to this work given that at least one planet in the vast majority 
of currently observed 2-planet systems will likely be affected by tides during post-MS evolution. 
Therefore understanding the fate of multi-planet systems crucially depends on effects such as 
the one we have modelled here. 



5 Conclusion 

We have explored properties of the three-body problem with mass loss. An astrophysical appli- 
cation for this work is two-planet post-MS systems. We proposed a mechanism, which combines 
the effects of stellar mass loss and mutual planetary gravitational interactions, that can desta- 
bilize a planetary system. Consequently, the lighter planet escapes and becomes an orphan 
planet. 

In particular, the coupling between mass loss and mutual interactions between the planets 
is not well-understood and was previously largely unexplored. We have identified the regions of 

2 See the Extrasolar Planet Encyclopedia at |http:/ /exoplanet.eu/| 
3 See the Exoplanet Data Explorer at http://exoplanets.org/ 
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Figure 17: The semimajor axis and eccentricity evolution of a system with mo = 4M , mi = 
0.6Mj and mi = 3Mj. The star loses mass linearly with rate a = 5 x 10~ 6 and (3 = 76%. 
The initial conditions are a\ = 10.0 and 02 = 18.0, e% = 0.1, e2 = 0.2, wi = 0°, il>2 = 180°, 
Mi = M 2 = 0°. 



phase space where the coupling is strong, causing initially stable systems to remain stable during 
the mass loss but later exhibiting escape. The stellar mass loss may transfer the system from 
its initially stable region to a wide chaotic sea in the phase space, which is associated with the 
final stellar mass. In such a chaotic region escape of the lighter planet always occurs. Although 
we find this evolutionary sequence to be robust to small changes in initial conditions, we find 
that the escape time is sensitive to these changes. However, it seems that chaotic evolution, 
which is detected by the maximal LCN during an integration interval up to ~ 100 ky, results in 
most cases to destabilization of the system in less than 100 My. The general phase space maps 
presented can be used to characterize a wide variety of real systems. 

In the numerical simulations of our paper we used a constant rate of mass loss. However, 
we have obtained qualitatively similar results for an exponential rate of mass loss: Fig. 4 
demonstrates excellent agreement between the constant mass loss model and the exponential 
model. Also, we note that mass loss along the Asymptotic Giant Branch closely resembles either 
constant or exponential mass loss. 

We find that the upper-mass end of potential planetary systems that will experience a white 
dwarf phase (M ~ 4 — 8M Q ) is likely to exhibit the chaotic evolutionary behaviour seen here. 
Just as importantly, lower-mass planetary systems, such as the currently-observed bound multi- 
planet exoplanet population, will experience adiabatic planet-planet interactions during post-MS 
mass loss. Claiming that the adiabatic approximation holds in these cases for the three-body 
problem is useful, and will simplify the analytics for future post-MS three-body studies. 
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